# Data manipulation
# ==============================================================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px

pd.set_option("display.max_columns", None)

# Data visualisation
# ==============================================================================
import geopandas as gpd
import osmnx as ox
import networkx as nx
import folium


# Load the GeoPackage file and check its structure
# ==============================================================================
gdf = gpd.read_file('distributioncenters2021NL_repo.gpkg')
display(gdf.shape)
display(gdf.columns)
display(gdf.head())

/opt/anaconda3/lib/python3.12/site-packages/pyogrio/geopandas.py:265: UserWarning: More than one layer found in 'distributioncenters2021NL_repo.gpkg': 'buildings2021NL' (default), 'estates2021NL', 'corridor2021NL'. Specify layer parameter to avoid this warning.
  result = read_func(

(26951, 28)

Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
       'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
       'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
       'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
       'Corop', 'Corridor', 'Current', 'geometry'],
      dtype='object')

# Filter the 'gdf' dataset to include 'Transport and Logistic services' and the city of Rotterdam, since it will be the focus of this analysis
# ==============================================================================
rotterdam = gdf[
(gdf['Function'] == 'Transport and logistic services') &
(gdf['Municip'] == 'Rotterdam')
]

# Check new dataframe structure
# ==============================================================================
display(rotterdam.shape)
display(rotterdam.head())
display(rotterdam.tail())

# Define the geographical location to download road network
# ==============================================================================
place_name = 'Rotterdam, Netherlands'

# Download the road network
# ==============================================================================
G = ox.graph_from_place(place_name, network_type= 'drive')

# Plot the street network
# ==============================================================================
ox.plot_graph(G)

(323, 28)

(<Figure size 800x800 with 1 Axes>, <Axes: >)

from scipy.spatial import cKDTree

# Convert the graph to GeoDataFrames
# ==============================================================================
nodes, edges = ox.graph_to_gdfs(G)

# Reproject Transport and logistics buildings to match network CRS
# ==============================================================================
rotterdam = rotterdam.to_crs(nodes.crs)

# Build KDTree from network node coordinates
# ==============================================================================
node_coords = np.array([(point.y, point.x) for point in nodes.geometry])
tree = cKDTree(node_coords)

# Get coords of logistics points (using centroids of polygons)
# ==============================================================================
building_coords = np.array([(geom.centroid.y, geom.centroid.x) for geom in rotterdam.geometry])

# Find nearest node index for each index
# ==============================================================================
distances, indices = tree.query(building_coords, k=1)

# Add nearest node ID to the 'rotterdam' dataframe
# ==============================================================================
rotterdam['nearest_node'] = nodes.iloc[indices].index

display(rotterdam.head(3))
display(rotterdam.columns)

Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
       'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
       'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
       'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
       'Corop', 'Corridor', 'Current', 'geometry', 'nearest_node'],
      dtype='object')

# Define warehouse origin and end point by index (select two random warehouses)
# ==============================================================================
origin_node = rotterdam.iloc[0]['nearest_node']
destination_node = rotterdam.iloc[20]['nearest_node']

# Convert node IDs to lat/lon coordinates
# ==============================================================================
origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
destination_coords = (G.nodes[destination_node]['y'], G.nodes[destination_node]['x'])

# Find the shortest path and the shortest path length using NetworkX (apply Dijkstra's algorithm)
# ==============================================================================
shortest_path = nx.shortest_path(G, origin_node, destination_node, weight= 'length')
path_length = nx.shortest_path_length(G, origin_node, destination_node, weight='length')

print(f" Shortest Path: {path_length:.2f} methers")

# Visualise using Folium 
# ==============================================================================
m = folium.Map(location=origin_coords, zoom_start=14)

folium.Marker(
    location=origin_coords, 
    icon=folium.Icon(color='blue', icon='info_sign')
).add_to(m) # Marker for origin point

folium.Marker(
    location=destination_coords, 
    icon=folium.Icon(color='green', icon='info_sign')
).add_to(m) # Marker for destination point

route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in shortest_path]
folium.PolyLine(route_coords, color='red', weight=4, opacity=0.8).add_to(m) # Add the route line

# Display the map
# ==============================================================================
m

 Shortest Path: 17558.97 methers

from shapely.geometry import Point

# Get bounding box coordinates of all logistics centers in Rotterdam
# (used to generate random delivery points within the city limits)
# ==============================================================================
minx, miny, maxx, maxy = rotterdam.total_bounds

# Generate 10 random points within the bounding box
# ==============================================================================
num_points = 10
x_coords = np.random.uniform(minx, maxx, num_points)
y_coords = np.random.uniform(miny, maxy, num_points)
delivery_points = [Point(x, y) for x, y in zip(x_coords, y_coords)]

# Create a GeoDataFrame for delivery points
# ==============================================================================
deliveries = gpd.GeoDataFrame(geometry=delivery_points, crs=rotterdam.crs)

# Initialize the map centered around Rotterdam
# ==============================================================================
m = folium.Map(location=[(miny + maxy) / 2, (minx + maxx) / 2], zoom_start=12)

# Build a KDTree of logistics center node coordinates
# ==============================================================================
center_nodes = rotterdam['nearest_node'].values
center_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in center_nodes])
center_tree = cKDTree(center_coords)

# Assign nearest road node to each delivery point
# ==============================================================================
delivery_coords = np.array([(point.x, point.y) for point in deliveries.geometry])
node_ids = list(G.nodes)
node_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in node_ids])
node_tree = cKDTree(node_coords)

# Match each delivery point to its nearest road node
# ==============================================================================
_, node_indices = node_tree.query(delivery_coords, k=1)
deliveries['delivery_node'] = [node_ids[i] for i in node_indices]

# Match each delivery to the closest logistics center based on network distance
# ==============================================================================
delivery_routes = []
for delivery_node in deliveries['delivery_node']:
    # Use Dijkstra to all centers
    lengths = nx.single_source_dijkstra_path_length(G, delivery_node, weight='length')
    
    # Pick closest center
    min_dist = float('inf')
    closest_center = None
    for center_node in center_nodes:
        if center_node in lengths and lengths[center_node] < min_dist:
            min_dist = lengths[center_node]
            closest_center = center_node
    
    delivery_routes.append({
        'delivery_node': delivery_node,
        'center_node': closest_center,
        'distance': min_dist
    })

# Add logistics center to the map
# ==============================================================================
for idx, row in rotterdam.iterrows():
    center_coords = (row.geometry.centroid.y, row.geometry.centroid.x)
    folium.Marker(location=center_coords, icon=folium.Icon(color='blue', icon='industry')).add_to(m)

# Add delivery points and routes to the map
# ==============================================================================
for idx, row in deliveries.iterrows():
    delivery_coords = (row.geometry.y, row.geometry.x)
    folium.Marker(location=delivery_coords, icon=folium.Icon(color='green', icon='home')).add_to(m)
    
    # Find the corresponding route
    route_info = delivery_routes[idx]
    path = nx.shortest_path(G, route_info['delivery_node'], route_info['center_node'], weight='length')
    route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in path]
    folium.PolyLine(route_coords, color='red', weight=2.5, opacity=0.8).add_to(m)


# Display the map
# ==============================================================================
m

# Choose a few random logistics centers (Origin Nodes)
# ==============================================================================
origin_nodes = rotterdam.iloc[[5, 24, 70]]['nearest_node'].values

colors = ['red', 'green', 'blue']
service_dists = [10000]  # or multiple if you want [5000, 10000]

m = folium.Map(location=origin_coords, zoom_start=11)

# Loop over each logistics center to generate and visualize its 10 km service area
# ==================================================================================
for idx, origin_node in enumerate(origin_nodes):
    for dist in service_dists:
        subgraph = nx.ego_graph(G, origin_node, radius=dist, distance='length')
        edges = ox.graph_to_gdfs(subgraph, nodes=False, edges=True)
        
        for _, row in edges.iterrows():
            line = folium.PolyLine(
                locations=[(y, x) for x, y in row['geometry'].coords],
                color=colors[idx],
                weight=2.5,
                opacity=0.5
            )
            line.add_to(m)

    # Add a marker for the center
    # ==============================================================================
    origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
    folium.Marker(location=origin_coords, icon=folium.Icon(color=colors[idx], icon='industry')).add_to(m)

m

import osmnx as ox
from pulp import *

# Download the Road Network
# ==============================================================================
city = "Rotterdam, Netherlands"
G = ox.graph_from_place(city, network_type="drive")
nodes, edges = ox.graph_to_gdfs(G)

np.random.seed(42)

# Sample demand points (e.g., customers)
# ==============================================================================
demand_points = nodes.sample(10).reset_index()
demand_points['id'] = ['D' + str(i) for i in range(len(demand_points))]

# Sample facility candidates (e.g., warehouses)
# ==============================================================================
facility_points = nodes.sample(5).reset_index()
facility_points['id'] = ['F' + str(i) for i in range(len(facility_points))]

# Calculate Distance Matrix
# ==============================================================================
demand_coords = demand_points[['y', 'x']].to_numpy()
facility_coords = facility_points[['y', 'x']].to_numpy()

distance_matrix = np.round(
    np.linalg.norm(demand_coords[:, None] - facility_coords, axis=2), 2
)

# Build and Solve the Optimization Model
# ==============================================================================
model = LpProblem("Facility_Location_Problem", LpMinimize)

# Decision variable: assign[d][f] = 1 if demand d assigned to facility f
# ==============================================================================
assign = LpVariable.dicts("Assign", 
                          [(d, f) for d in range(len(demand_points)) for f in range(len(facility_points))],
                          cat='Binary')

# Open facility variable: 1 if facility f is open
# ==============================================================================
open_facility = LpVariable.dicts("Open", 
                                 list(range(len(facility_points))), 
                                 cat='Binary')

# Objective: Minimize total distance
# ==============================================================================
model += lpSum(assign[d, f] * distance_matrix[d][f]
               for d in range(len(demand_points))
               for f in range(len(facility_points)))

# Constraint: Each demand must be assigned to one facility
# ==============================================================================
for d in range(len(demand_points)):
    model += lpSum(assign[d, f] for f in range(len(facility_points))) == 1

# Constraint: Only assign to open facilities
# ==============================================================================
for d in range(len(demand_points)):
    for f in range(len(facility_points)):
        model += assign[d, f] <= open_facility[f]

# Limit number of open facilities
# ==============================================================================
model += lpSum(open_facility[f] for f in range(len(facility_points))) <= 2

# Solve it
# ==============================================================================
model.solve()

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /opt/anaconda3/lib/python3.12/site-packages/pulp/apis/../solverdir/cbc/osx/i64/cbc /var/folders/p4/km2q7qsd1cl5b1h60220jf8w0000gn/T/f29e0b4fc3e04cc7a473591e24688089-pulp.mps -timeMode elapsed -branch -printingOptions all -solution /var/folders/p4/km2q7qsd1cl5b1h60220jf8w0000gn/T/f29e0b4fc3e04cc7a473591e24688089-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 66 COLUMNS
At line 382 RHS
At line 444 BOUNDS
At line 500 ENDATA
Problem MODEL has 61 rows, 55 columns and 155 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 0.36 - 0.00 seconds
Cgl0004I processed model has 61 rows, 55 columns (55 integer (55 of which binary)) and 155 elements
Cutoff increment increased from 1e-05 to 0.00999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 0.36
Cbc0038I Before mini branch and bound, 55 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (0.01 seconds)
Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 0.36 - took 0.00 seconds
Cbc0012I Integer solution of 0.36 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds)
Cbc0001I Search completed - best objective 0.36, took 0 iterations and 0 nodes (0.01 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 0.36 to 0.36
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                0.36000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.02

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.02

1

# Extract coordinates
# ==============================================================================
demand_coords = demand_points[['y', 'x']].to_numpy()
facility_coords = facility_points[['y', 'x']].to_numpy()

# Identify which facilities are open
# ==============================================================================
open_facilities = [f for f in range(len(facility_points)) if open_facility[f].varValue == 1]

# Create a plot
# ==============================================================================
fig, ax = plt.subplots(figsize=(10, 10))

# Plot the road network as background
# ==============================================================================
edges.plot(ax=ax, linewidth=0.2, color='gray', alpha=0.4)

# Plot demand points
# ==============================================================================
plt.scatter(demand_coords[:, 1], demand_coords[:, 0], c='red', label='Demand Points')

# Plot opened facilities
# ==============================================================================
for f in open_facilities:
    plt.scatter(facility_coords[f][1], facility_coords[f][0], c='blue', s=100, label='Opened Facility' if f == open_facilities[0] else "")

# Draw assignment lines
# ==============================================================================
for d in range(len(demand_points)):
    for f in open_facilities:
        if assign[d, f].varValue == 1:
            x = [demand_coords[d][1], facility_coords[f][1]]
            y = [demand_coords[d][0], facility_coords[f][0]]
            plt.plot(x, y, 'g--', linewidth=1)

# Add legend and title
# ==============================================================================
plt.legend()
plt.title("Facility Location Optimization: Opened Facilities and Assignments")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.grid(True)
plt.show()

	osm_id	BAG_id	Constr_yea	Year_cat	Footprint	Size_cat	Function	Function_n	Job_dens	Job_dens_n	RIN	Year	Name	Locality	Gross_ha	Net_ha	Max_cat	Rail_cat	Water_cat	Muni_nr	Municip	Prov_nr	Province	Corop_nr	Corop	Current	geometry
0	273364927	0441100000002816	1972	1970s	778.057432	S	Trade, Import, Export	1	medium	3	1741	1971.0	Witte Paal	Schagen	76.0000	43.0000	4	A	E	441	Schagen	27	Noord-Holland	18	Kop Van Noord-holland	0	MULTIPOLYGON Z (((115780.809 532845.953 0, 115...
1	279126121	0971100000022139	1978	1970s	734.322479	S	Trade, Import, Export	1	medium	3	203678	1960.0	Kerensheide	Stein	19.8223	16.9625	3	D	E	971	Stein	31	Limburg	39	Zuid-limburg	1	MULTIPOLYGON Z (((182791.839 331932.418 0, 182...
2	615938117	1931100000000715	2018	2010s and later	836.482827	S	Trade, Import, Export	1	low	2	171034	2012.0	Zevender	Schoonhoven	15.0000	10.7000	4	A	E	1931	Krimpenerwaard	28	Zuid-Holland	28	Oost-zuid-holland	1	MULTIPOLYGON Z (((118809.492 440320.413 0, 118...
3	284074624	0302100000119189	2013	2010s and later	2048.445693	S	XXL retail center	3	low	2	1203	1943.0	Kom Nunspeet	Nunspeet	8.7132	8.3971	5	A	E	302	Nunspeet	25	Gelderland	13	Veluwe	1	MULTIPOLYGON Z (((181723.773 487021.56 0, 1817...
4	280745712	0050100000352829	1999	1990s	1610.136749	S	Trade, Import, Export	1	very high	5	213	1975.0	Trekkersveld I en II	Zeewolde	160.1670	93.9994	5	A	A	50	Zeewolde	24	Flevoland	40	Flevoland	1	MULTIPOLYGON Z (((163088.239 484983.775 0, 163...

	osm_id	BAG_id	Constr_yea	Year_cat	Footprint	Size_cat	Function	Function_n	Job_dens	Job_dens_n	RIN	Year	Name	Locality	Gross_ha	Net_ha	Max_cat	Rail_cat	Water_cat	Muni_nr	Municip	Prov_nr	Province	Corop_nr	Corop	Corridor	Current	geometry
33	271881098	0599100000653535	1999	1990s	784.190901	S	Transport and logistic services	2	very high	5	170371	2010.0	Europoort	Rotterdam	1700.0	1369.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((72649.646 436919.835 0, 7266...
75	273554871	0599100000643345	1990	1990s	2362.159063	S	Transport and logistic services	2	very high	5	170143	1964.0	Prins Alexander	Rotterdam-noord	35.0	29.0	3	A	E	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((97149.757 440577.129 0, 9718...
181	554042637	None	1994	1990s	4874.941239	M	Transport and logistic services	2	no data	99	170380	1973.0	Eemhaven	Rotterdam	358.0	292.0	4	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((88339.665 433065.91 0, 88329...
228	273326689	0599100010036685	1988	1980s	3352.465527	M	Transport and logistic services	2	no data	99	170373	1970.0	Botlek Europoort-oost	Rotterdam	1200.0	1069.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((77162.415 434307.927 0, 7718...
229	273326686	0599100000754059	2001	2000s	10457.050191	L	Transport and logistic services	2	low	2	170376	1990.0	Distripark Botlek	Rotterdam	120.0	108.0	4	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((76763.925 432037.599 0, 7677...

	osm_id	BAG_id	Constr_yea	Year_cat	Footprint	Size_cat	Function	Function_n	Job_dens	Job_dens_n	RIN	Year	Name	Locality	Gross_ha	Net_ha	Max_cat	Rail_cat	Water_cat	Muni_nr	Municip	Prov_nr	Province	Corop_nr	Corop	Corridor	Current	geometry
26041	273554710	0599100000617012	1976	1970s	1782.039078	S	Transport and logistic services	2	very low	1	170143	1964.0	Prins Alexander	Rotterdam-noord	35.0	29.0	3	A	E	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((96511.132 440308.944 0, 9654...
26279	863262550	0599100110003461	2020	2010s and later	33920.430577	XL	Transport and logistic services	2	no data	99	170341	2008.0	Maasvlakte	Rotterdam	2728.0	2297.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((58428.859 439331.876 0, 5851...
26280	863262551	0599100110004080	2020	2010s and later	20116.635052	XL	Transport and logistic services	2	no data	99	170341	2008.0	Maasvlakte	Rotterdam	2728.0	2297.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((58522.069 439072.086 0, 5864...
26445	271523822	0599100000659860	2000	2000s	5113.122114	M	Transport and logistic services	2	very high	5	170164	1971.0	Hoogvliet	Hoogvliet	28.0	19.8	3	A	E	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((85042.696 432049.426 0, 8506...
26896	271871679	0599100000047089	1987	1980s	803.993757	S	Transport and logistic services	2	very high	5	170371	2010.0	Europoort	Rotterdam	1700.0	1369.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((71684.971 437914.726 0, 7170...

	osm_id	BAG_id	Constr_yea	Year_cat	Footprint	Size_cat	Function	Function_n	Job_dens	Job_dens_n	RIN	Year	Name	Locality	Gross_ha	Net_ha	Max_cat	Rail_cat	Water_cat	Muni_nr	Municip	Prov_nr	Province	Corop_nr	Corop	Corridor	Current	geometry	nearest_node
33	271881098	0599100000653535	1999	1990s	784.190901	S	Transport and logistic services	2	very high	5	170371	2010.0	Europoort	Rotterdam	1700.0	1369.0	5	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((4.19026 51.91467 0, 4.19049 ...	44519439
75	273554871	0599100000643345	1990	1990s	2362.159063	S	Transport and logistic services	2	very high	5	170143	1964.0	Prins Alexander	Rotterdam-noord	35.0	29.0	3	A	E	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((4.5457 51.95062 0, 4.54624 5...	44550916
181	554042637	None	1994	1990s	4874.941239	M	Transport and logistic services	2	no data	99	170380	1973.0	Eemhaven	Rotterdam	358.0	292.0	4	B	A	599	Rotterdam	28	Zuid-Holland	29	Groot-rijnmond	1	1	MULTIPOLYGON Z (((4.41902 51.88213 0, 4.41887 ...	1582380654

Optimizing Urban Logistics with Network Analysis in Rotterdam¶

Shortest Possible Route¶

Closest Facility Analysis¶

Service Area Analysis¶

Facility Location Optimization¶

Conclusions¶